#include <iostream>
#include <Eigen/Dense>

int test_main2();
int test_main3();
int test_main4();

void ComputeCoeffs(Eigen::VectorXd& x, Eigen::VectorXd& y, Eigen::VectorXd& coeffs)
{
    // 构建范德蒙矩阵
    Eigen::MatrixXd A(x.size(), 3);
    for (int i = 0; i < x.size(); i++) {
        A(i, 0) = x[i] * x[i]; // x^2
        A(i, 1) = x[i];       // x
        A(i, 2) = 1;          // 1
    }

    // 最小二乘求解
    coeffs = A.householderQr().solve(y);
    std::cout << "拟合系数: [a2 a1 a0] = " << coeffs.transpose() << std::endl;
}

Eigen::VectorXd polyfit(const Eigen::VectorXd& x, const Eigen::VectorXd& y, int degree) {
    Eigen::MatrixXd A(x.size(), degree + 1);
    for (int i = 0; i < x.size(); i++) {
        for (int j = 0; j <= degree; j++) {
            A(i, j) = pow(x(i), j);  // 构建范德蒙矩阵
        }
    }
    return A.householderQr().solve(y);  // QR分解求解
}

Eigen::VectorXd ComputeFitValue(double value, Eigen::VectorXd &coeffs)
{
    Eigen::VectorXd z(1);

    Eigen::MatrixXd B(z.size(), 3);

    z << value;
    for (int i = 0; i < z.size(); i++) {
        B(i, 0) = z[i] * z[i]; // x^2
        B(i, 1) = z[i];       // x
        B(i, 2) = 1;          // 1
    }
    // 计算拟合值
    Eigen::VectorXd z_fit = B * coeffs;
    return z_fit;
}

int main() 
{
    test_main4();
    test_main3();
    test_main2();

    Eigen::VectorXd x(5);
    Eigen::VectorXd y(5);
    x << 1, 2, 3, 4, 5;
    y << 1.1, 3.9, 8.2, 15.3, 24.5;
    //y << 3, 5, 8, 9, 12;

    // 构建范德蒙矩阵
    Eigen::MatrixXd A(x.size(), 3);
    for (int i = 0; i < x.size(); i++) {
        A(i, 0) = x[i] * x[i]; // x^2
        A(i, 1) = x[i];       // x
        A(i, 2) = 1;          // 1
    }

    Eigen::VectorXd coeffsy;
    ComputeCoeffs(x, y, coeffsy);

    // 最小二乘求解
    Eigen::VectorXd coeffs = A.householderQr().solve(y);
    std::cout << "拟合系数: [a2 a1 a0] = " << coeffs.transpose() << std::endl;
    std::cout << "The coefficients vector is: \n" << std::endl;
    for (int32_t i = 0; i < coeffs.size(); i++)
    {
        std::cout << "coeffs" << i << ": " << coeffs(i) << std::endl;
    }

    // 计算拟合值
    Eigen::VectorXd y_fit = A * coeffs;
    std::cout << "拟合结果:\n" << y_fit << std::endl;

    //Eigen::VectorXd z(1);

    //Eigen::MatrixXd B(z.size(), 3);
    //z << 1.5;
    //for (int i = 0; i < z.size(); i++) {
    //    B(i, 0) = z[i] * z[i]; // x^2
    //    B(i, 1) = z[i];       // x
    //    B(i, 2) = 1;          // 1
    //}
    //// 计算拟合值
    //Eigen::VectorXd z_fit = B * coeffs;
    //std::cout << "拟合结果:\n" << z_fit << std::endl;

    Eigen::VectorXd aa_fit = ComputeFitValue(1.5, coeffs);
    std::cout << "拟合结果:\n" << aa_fit[0] << std::endl;

    Eigen::VectorXd coeffs2 = polyfit(x, y, 2);
    std::cout << "拟合系数: " << coeffs2.transpose() << std::endl;

    return 0;
}